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Abstract 


In  this  paper  we  present  a  routine  that  exploits  the  power  of  seismic  arrays  and  cepstral 
techniques  to  estimate  the  depth  of  an  event  directly  from  the  observed  seismograms.  A 
discussion  of  the  pertinent  geophysical  assumptions,  cepstral  processing  algorithm,  stable  peak 
identification  via  “cepstrograms,”  false  alarm  reduction  methodology,  and  our  array-based  depth 
estimation  routine  is  presented.  An  analysis  of  several  shallow  events  is  performed  and 
compared  to  results  produced  by  a  standard  location  algorithm,  waveform  forward  modeling,  and 
previously  published  solutions. 
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Automated  Source  Depth  Estimation  Using  Array  Processing  Techniques 


1.0  Introduction 

Source  depth  estimation  is  a  key  process  in  the  discrimination  of  earthquakes  and  explosions. 
The  lack  of  observable  depth  phases  does  not  necessarily  mean  the  event  occurred  at  or  near  the 
surface.  Shallow  events  can  have  closely  spaced  depth  phases  that  are  indistinguishable  even  by 
seasoned  human  analysts.  Moreover,  the  onset  of  smaller  events  observed  at  regional  distances 
is  often  complicated  by  the  arrival  of  multiple  phases  in  rapid  succession,  which  makes  the 
identification  of  depth  phases  even  more  problematic.  Source  parameters  for  such  events  can  be 
derived  using  moment  tensor  inversion  or  forward  modeling  techniques,  which  are  difficult  to 
apply  to  events  less  than  mb  5.5  and  shallower  than  15  km,  and  depend  on  the  availability  and 
accuracy  of  geophysical  models.  These  limitations  are  not  practical  for  real-time  discrimination 
of  earthquakes  and  explosions. 

If  depth  phases  with  sufficient  signal-to-noise  ratio  (SNR)  reside  in  an  observation,  they  will 
produce  a  spectral  scalloping  pattern  with  a  period  equal  to  the  time  delay  between  signals.  This 
spectral  phenomenon  can  be  detected  using  cepstral  processing,  which  has  been  used  in  a 
number  of  studies  over  the  last  45  years  with  limited  success  [Bogart  et  al.,  1963;  Ulrych,  1971; 
Kemerait  and  Childers,  1972;  Ulrych  et  al.,  1972;  Tribolet,  1978;  Kemerait  and  Sutton,  1982; 
Marenco  and  Madisetti,  1997;  Shumway  et  al.,  1998;  Bonner  et  al.,  2002;  and  Reiter,  2005]. 
These  studies,  however,  did  not  exploit  the  power  of  seismic  arrays  to  determine  the  ray 
parameter  of  the  arriving  phase.  The  ray  parameter,  an  assumed  wave  speed,  and  simple  vector 
decomposition  can  be  used  to  determine  the  vertical  phase  velocity  and  wavefront  angle  of 
incidence.  If  reciprocity  between  the  source  and  receiver  holds,  the  angle  of  incidence  and  take¬ 
off  angle  are  the  same  and  can  be  combined  with  the  depth  phase  delay  time  to  calculate  a  source 
depth  directly  from  the  observed  seismograms.  Unlike  moment  tensor  inversion  or  waveform 
forward  modeling,  this  methodology  neither  requires  detailed  geophysical  models  nor  is 
restricted  to  large  events  or  a  minimum  depth. 

Our  routine  employs  a  multi-stage  detection  scheme  that  reduces  the  high  false  alarm  rate 
inherent  to  cepstral  analysis.  First,  a  site-specific,  adaptive,  cepstral  amplitude  or  gamnitude 
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threshold,  recalling  the  terminology  coined  by  Bogart  et  al.,  1963,  is  derived  using  pre-signal 
noise  to  identify  statistically  significant  peaks.  Knowing  that  cepstra  are  highly  unstable  and 
change  significantly  with  minor  changes  to  processing  parameters,  we  developed  an  iterative 
technique  to  search  for  stable  detections  over  a  series  of  increasing  time  windows.  The  resulting 
“cepstrograms”  accentuate  stable  features  in  the  cepstral  domain  to  assist  the  algorithm  in 
selecting  only  signal-induced  peaks.  Finally,  a  binary  stacking  module  checks  for  consistent 
detections  across  the  observing  network. 

2.0  Theory 

Figure  la  shows  the  ray  path  geometry  between  a  shallow  source  and  receiver.  Notice  that 
upward  traveling  rays  reflect  off  the  free  surface,  travel  along  a  path  similar  to  the  primary  phase, 
and  arrive  at  the  receiver  with  the  same  angle  of  incidence  and  apparent  velocity.  This  idealized 
illustration  depicts  the  geophysical  assumptions  our  algorithm  relies  on,  which  are  as  follows: 

•  Source  Mechanism:  Cepstral  analysis  relies  on  the  assumption  that  the  source 
mechanism  can  be  modeled  as  a  point  source.  Large  magnitude  earthquakes  often  have 
time  varying  rupture  processes  that  violate  this  assumption.  As  a  result,  we  limit 
ourselves  to  analyzing  events  with  bodywave  magnitudes  less  than  6.0. 

•  Phase  Speed:  Since  we  are  interested  in  discriminating  between  earthquakes  and 
explosions,  we  assume  a  shallow  source  depth  (d  <  20  km).  This  means  that  the  speed  of 
the  incident  P-wave  is  approximately  5.8  km/sec  for  continental  crust  events  [Kennett, 
1991]. 

•  Angle  of  Arrival:  If  reciprocity  holds,  the  incidence  angle  of  the  primary  arrival  is  equal 
to  the  take-off  angle  at  the  source.  This  assumption  allows  for  the  derivation  of  the  take¬ 
off  angle  using  horizontal  apparent  velocity  measurements  and  the  previously  assumed 
phase  speed. 
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Receiver 


(a) 


(b) 


Figure  1.  (a)  Propagation  paths  for  primary  arrival  and  associated  depth  phases  traveling  through  a  flat, 
discretized  earth,  (b)  Ideal  ray  path  geometry  for  a  seismic  point  source. 


2.1  Depth  Estimation 

The  ray  transmission  and  reflection  geometry  generated  by  a  seismic  point  source  is  shown  in  Figure 
lb.  The  illustration  shows  depth,  d  (km),  is  a  function  of  the  delay  time  between  the  primary  arrival  and 
its  associated  depth-phase,  x  (sec),  the  ray  take-off  angle,  9  (deg),  and  the  P  wave  speed,  a  (km/s) 

d  =  (T*a)  r  — — 1 

L2cos6^J  Qx 


The  value  of  x  is  supplied  via  cepstral  processing  (section  2.2.1)  and  the  ray  take-off  angle  is  computed 
using  the  phase  velocity's  horizontal  and  vertical  components.  The  horizontal  phase  velocity  or  apparent 
velocity,  c  (s/km),  of  a  planar  wavefront  traveling  across  a  seismic  array  is  often  measured  using 
frequency-wavenumber  analysis  [Kvaerna,  1989].  The  apparent  velocity  measurement  of  the  incident 
wavefront,  and  an  assumed  speed  of  the  P-wave,  allows  us  to  calculate  the  ray's  vertical  velocity 
component,  r|  (s/km),  and  take-off  angle  using 


respectively. 


1 


0  =  tan 


(2) 

(3) 
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Since  the  take-off  angle  and  apparent  velocity  vary  as  a  function  of  distance,  due  to  varying  ray  path 
geometries,  they  are  calculated  for  each  station  in  the  network.  Site-specific  values  for  x,  c,  and  0  are 
computed  and  substituted  into  (1)  [Junek  et  ah,  2006;  Junek  et  ah,  2007].  This  results  in  a  suite  of  depth 
hypothesis  that  already  account  for  move-out  between  stations. 

2.2  Signal  Processing  Algorithm 

Our  signal  processing  algorithm  (Figure  2)  consists  of  three  main  components.  First,  a  cepstral 
processing  component  (cyan)  determines  the  time  delays  between  direct  P  and  the  reflected  phase.  Next, 
a  depth  estimation  component  (yellow)  combines  delay  times  with  phase  wavespeeds  to  compute  a  depth. 
Finally,  a  false  alarm  reduction  component  (green)  identifies  statistically  significant  cepstra  and 
consistent  depth  estimates  across  the  network. 


Figure  2.  Signal  processing  algorithm  consists  of  a  cepstral  processing  component  (cyan),  depth 
estimation  component  (yellow),  and  a  false  alarm  reduction  component  (green). 


4 


Automated  Source  Depth  Estimation 
Using  Array  Processing  Techniques 

2.2.1  Cepstral  Processing 

Our  cepstral  processing  function  combines  the  methodologies  of  [Kemerait,  1972;  Shumway  et  ah, 
1998;  and  Bonner  et  ah,  2002].  Event  observations  and  pre-signal  noise  segments  from  each  element  of 
an  array  are  passed  into  the  algorithm  and  processed  separately  using  the  same  parameters  to  ensure  the 
results  are  comparable.  Mean  signal  and  pre-signal  noise  cepstra  are  created  from  the  individual  results  to 
enhance  common  peaks. 

2.2.2  False  Alarm  Reduction 

The  false  alarm  reduction  routine  consists  of  four  primary  components:  gamnitude  threshold 
computation,  detection  processing,  application  of  a  cepstral  stability  requirement,  and  a  network 
consistency  check.  Each  of  these  techniques  is  used  to  reduce  the  high  false  alarm  rate  inherent  to 
cepstral  analysis. 

A  gamnitude  threshold  derived  from  site-specific,  pre-signal  noise  cepstra  is  used  to  select  candidate 
peaks  for  the  depth  estimation  algorithm.  The  threshold  is  defined  as  the  99^^  percentile  of  the  gamnitude 
distribution  of  the  pre-signal  noise  cepstra  for  a  sampling  window  equal  to  the  time -domain  sampling 
window.  This  is  repeated  for  each  station  in  the  network  to  derive  real-time,  site-specific  gamnitude 
thresholds  that  are  based  on  the  current  noise  conditions  at  each  site.  This  prevents  hourly,  daily,  or 
seasonal  noise  fluctuations  from  increasing  the  false  alarm  rate. 

Cepstral  processing  is  performed  for  each  station  using  a  series  of  increasing  sampling  window 
lengths  to  identify  stable  peaks.  As  the  sampling  window  length  grows  and  captures  larger  sections  of  the 
depth  phase,  the  intensity  of  the  points  in  the  ‘‘cepstrograms”  grows,  peaks,  and  fades  as  more  noise  is 
acquired.  A  stability  parameter,  k,  is  used  to  define  the  number  of  consecutive  threshold  crossing  cepstra 
that  are  required  to  declare  candidate  depth  phase  detections.  The  value  of  k  is  typically  set  between  15% 
and  25%  of  the  total  number  of  sampling  windows.  Results  existing  for  less  than  this  value  are  not 
considered  a  candidate  depth  phase.  Figure  3  shows  conceptual  cepstrograms  for  both  pre-signal  noise 
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and  the  observed  seismogram.  This  proeedure  is  repeated  for  eaeh  set  of  array  observations  until  a 
eolleetion  of  eepstrograms  are  generated. 
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Figure  3.  Coneeptual  eepstrograms  for  pre-signal  noise  and  observed  seismograms,  respeetively,  where 
the  Y-axis  is  delay  time,  X-axis  is  the  sampling  window  length,  random  points  in  the  top  and 
bottom  panels  are  transient  noise  spikes,  and  lines  are  stable  signals. 


Deteetion  proeessing  is  earried  out  on  a  station-by-station  basis.  All  threshold  erossing  eepstra,  for 
eaeh  station,  that  meet  the  stability  eriteria  are  treated  equally  to  avoid  the  possibility  of  a  missed 
deteetion.  Eaeh  eandidate  depth  phase  delay  time  is  then  passed  to  the  depth  estimator  (seetion  2.2.3). 

Network  eonsisteney  is  eheeked  by  a  binary  staeking  algorithm  and  takes  plaee  after  depth  extraetion 
to  eompensate  for  move-out  between  stations  [Murphy  et  al.,  1999;  Bonner  et  al.,  2002].  This 
methodology  allows  one  input  per  station  for  eaeh  depth  eell,  whose  width  is  a  user  defined  parameter,  n 
[Bonner  et  al.,  2002;  Murphy  et  al.  1999].  The  largest  peak  in  the  staek  identifies  the  measurement  that  is 
the  most  eonsistent  aeross  the  network  and  is  deelared  the  final  result. 
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2.2.3  Depth  Computation 

The  depth  estimation  module  requires  time  domain  data  for  the  primary  arrival  and  the  time  delay  of 
eaeh  threshold  erossing  eepstra  for  eaeh  station  being  eonsidered.  Time  domain  data  is  used  to  eompute  e 
of  the  ineident  wavefront,  whieh  is  used  to  eompute  p  and  9.  These  parameters  are  substituted  into 
equation  (1)  to  ealeulate  a  suite  depth  estimates  for  that  seismogram.  This  is  repeated  for  eaeh  array  in 
the  network,  where  the  resulting  depth  profiles  are  passed  to  the  network  eonsisteney  routine  (seetion 
2.2.2). 

3.0  Discussion 

Automated  depth  estimates  for  a  series  of  events  observed  at  regional  and  teleseismie  distanees  were 
generated  and  eompared  to  those  derived  by  a  standard  loeation  routine,  moment  tensor  inversion  and 
waveform  forward  modeling,  and  previously  published  results.  Five  randomly  seleeted  events  are  ehosen 
for  our  evaluation.  Filter  passbands  that  maximize  the  signal-to-noise  ratio  for  eaeh  station/event  pair 
were  seleeted  and  a  standard  set  of  proeessing  parameters  were  used  to  prevent  tuning  biases  in  the 
solutions. 

Models  for  events  1,  2,  4,  and  5  were  eomputed  using  the  Moment  Tensor  Inversion  Toolkit 
(MTINV)  [lehinose,  2006]  and  regional  data  aequired  from  IRIS  or  the  Japanese  Meteorologieal  Ageney 
(JMA).  Simple  three-layer  erustal  models  over  a  half  spaee  were  used  to  model  these  events,  where  the 
Western  United  States  model  was  used  for  events  2  and  5,  a  model  ereated  by  [lehinose,  2008]  was  used  for 
event  4,  and  a  modified  version  of  a  model-based  one  [lehinose  et  ah,  2005]  was  used  for  event  1.  Event 
3  was  modeled  using  refieetivity  software  employing  Kennetf  s  teehnique  of  solving  wave  propagation 
problems  in  laterally  homogenous  layers  [Randall  et  ah,  no  date].  A  186-layer  Earth  model  eonsisting  of  a 
two-layer  erust,  similar  to  one  used  by  [Antolik  and  Drenger,  2003],  and  an  upper  and  lower  mantle 
model  based  on  PREM  was  used  to  eompute  the  syntheties  [Randall,  2006]. 

Observed  waveforms,  network  eonfiguration,  and  automated  proeessing  results  for  event  3  are  shown 
in  Figures  4  and  5.  Cepstra  for  eaeh  station  were  generated  for  a  series  of  sampling  window  lengths 
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between  0.0  see  and  8.0  see  in  inereasing  inerements  of  0.05  see.  Resulting  eepstrograms  and  frequeney 
wavenumber  plots  are  shown  in  Figures  5a  and  b,  respeetively.  Notiee  there  are  numerous  features  in 
eaeh  eepstrogram.  The  before  mentioned  stability  parameter  sereened  out  the  transients  and  passed  only 
stable  features  to  the  depth  eomputation  module.  The  final  depth  estimate  is  shown  in  Figure  5e  and  is 
approximately  3  km. 


Figure  4.  (a)  Event  #3:  2003  Bhuj  aftershoek  loeation,  foeal  meehanism,  and  network  eonfiguration. 
(b)  Seismograms  observed  by  eaeh  station  and  separated  as  a  funetion  of  distanee. 
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A 


Figure  5.  (a)  Cepstrograms  for  BURAR,  BRTR,  FINES,  and  KSRS,  where  the  top  and  bottom  panels 
for  eaeh  station  represent  pre-signal  noise  and  observed  seismogram  eepstra,  respeetively. 
The  floor  of  eaeh  eepstrogram  is  set  to  the  site-speeifle  gamnitude  threshold,  where  eaeh 
visible  feature  is  a  threshold  erossing  eepstra.  Delay  times  between  0  and  8  seeonds  were 
eonsidered  in  our  analysis;  however,  only  the  0-  to  2-seeond  delay  time  range  is  shown  for  the 
purpose  of  elarity.  (b)  Frequeney  wavenumber  plots,  (c)  Network-based  depth  estimate 
eorresponds  to  the  largest  peak,  whieh  is  approximately  3  km. 
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Table  1  shows  a  comparison  of  depth  estimates  for  the  analyzed  events.  Values  listed  in  the  ‘‘Array- 
Based”  column  were  produced  by  our  routine,  free-depth  solutions  were  created  using  a  location 
algorithm  based  on  [Jordan  and  Sverdrup,  1981],  and  published  solutions  were  obtained  from  several 
organizations,  which  are  referenced  in  Table  1.  Our  results  are  in  good  agreement  with  the  published  and 
free-depth  solutions  and  correspond  particularly  well  to  the  modeled  results. 


T able  1 :  Comparison  of  Results 


Depth  (km) 

Event 

Origin  Time  (GMT) 

Latitude  (°) 

Longitude  (°) 

mb 

Published 

Array- 

Based 

Modeled 

Free- 

Depth 

05/06/2002,  08:12:14 

38.4° 

141.2° 

5.1 

40 

35 

35 

45 

^  ** 

2 

03/12/2005,07:36:10 

39.2° 

40.8° 

5.4 

16 

5 

10 

13 

3** 

08/05/2003,  08:04:05 

23.7° 

70.4° 

5.1 

15 

3 

2^^ 

16 

4 

01/14/2004, 16:58:48 

27.5° 

52.17° 

5.4 

12 

2 

6 

10 

5" 

05/01/2006,  00:39:26 

w— - - - : - : - : — ^ - 

42.4° 

69.2° 

4.6 

15 

21 

18 

17 

Japanese  Metrological  Agency 
**Harvard  CMT 
+  KNDC  Solution 

++  Modeled  using  Reflectivity  Method 


4.0  Summary 

Our  automated,  array-based  depth  estimation  routine  produced  results  that  are  in  good  agreement  with 
those  created  by  conventional  methods.  The  false  alarm  reduction  processes  increased  the  reliability  of 
the  algorithm  by  selecting  cepstra  that  were  greater  than  or  equal  to  the  99^^  percentile  of  the  pre-signal 
noise  gamnitude  distribution  and  exist  across  multiple  sites.  Our  adaptive  detection  threshold  was  derived 
from  the  current  noise  conditions  at  each  site,  which  prevented  daily  noise  fluctuations  from  producing 
false  alarms.  Applying  the  stability  parameter  resulted  in  the  selection  of  highly  robust  features  in  the 
cepstrogram  and  screened  transient  noise  features  that  would  have  produces  false  depth  estimates. 
Moreover,  the  network  consistency  check  reduced  the  possibility  of  anomalous  cepstral  peaks  producing 
false  alarms  by  requiring  a  result  to  exist  across  multiple  sites. 
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The  combination  of  cepstral  processing  and  frequency  wavenumber  analysis  resulted  in  a  fast  and 
simple  technique  that  can  be  executed  in  near  real-time.  Unlike  moment  tensor  inversion  or  waveform 
modeling,  our  routine  requires  neither  detailed  geophysical  models  nor  is  restricted  to  large  events. 
Analysis  of  a  small  group  of  events  showed  its  ability  to  estimate  the  depth  of  extremely  shallow  events 
and  its  potential  as  a  real-time  discrimination  tool  for  cases  where  depth  phases  are  not  perceptible. 
Future  work  will  focus  on  applying  this  technique  to  larger  data  sets  and  the  routine  analysis  real-time 
data. 
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